#ifndef _DTI_FIBER_TRACKING_CPU
#define _DTI_FIBER_TRACKING_CPU

#include "../DWMRICore/DTI.h"
#include "../DWMRICore/Fiber.h"
#include "../DWMRICore/Scalar.h"

enum FT_METHOD {
	FT_STT,					/* basic streamline method */
	FT_STT_MLS,				/* moving least square streamline method */
	FT_FACT,				/* FACT fiber tracking method */
	FT_TEND,				/* tensor deflection fiber tracking method */
	FT_TEND_MLS				/* tensor deflection & moving least square */
};

class CDTIFiberTrackingCPU
{
public:
	/* input data */
	CDTI *m_pDTI;
	CScalar *m_pMask;

	CScalar *m_pFA;

	/* data dimensions */
	uint m_width;
	uint m_height;
	uint m_depth;

	/* general fiber tracking parameters */
	float m_maxLength;
	float m_stepSize;
	float m_faThreshold;
	float m_angleThreshold;

	/* MLS fiber tracking parameters */
	float m_MLSSize;
	CPoint3F *m_pSamplePoints;
	float *m_pSamplePointsWeight;

	/* tensor deflection fiber tracking parameters */
	float m_tendF;
	float m_tendG;

public:
	CDTIFiberTrackingCPU(void);
	~CDTIFiberTrackingCPU(void);

public:
	int InitDTIFiberTracking(const char *DTIPathname, const char *MaskPathname, 
		const float stepSize, const float maxLength,
		const float faThreshold, const float angleThreshold);

	void SetMLSParameters(const float size) { m_MLSSize = size; }
	void SetTENDParameters(const float f, const float g) { m_tendF = f; m_tendG = g; }

public:
	void FiberTracking(FT_METHOD method, const char *fiberPathname, const uint scale);
	
	/* basic streamline fiber tracking algorithm, implemented in FiberTracking_STT.cpp */
private:
	void FiberTracking_STT(const char *fiberPathname, const uint scale);
	void GenerateSeeds_STT(Seed *seeds, const uint scale, const int z);
	void TrackingOneFiber_STT(Seed *seed, CFiber *fiber);
	int GetNextDirection_STT(CPoint3F *pos, CVector3F *in, CVector3F *out);

	/* tensor deflection fiber tracking algorithm, implemented in FiberTracking_TEND */
private:
	void FiberTracking_TEND(const char *fiberPathname, const uint scale);
	void GenerateSeeds_TEND(Seed *seeds, const uint scale, const int z);
	void TrackingOneFiber_TEND(Seed *seed, CFiber *fiber);
	int GetNextDirection_TEND(CPoint3F *pos, CVector3F *in, CVector3F *out);

	/* tensor deflection with MLS filter fiber tracking algorithm, implementd in FiberTracking_TENDMLS.cpp */
private:
	void FiberTracking_TENDMLS(const char *fiberPathname, const uint scale);
	void GenerateSeeds_TENDMLS(Seed *seeds, const uint scale, const int z);
	void TrackingOneFiber_TENDMLS(Seed *seed, CFiber *fiber);
	int GetNextDirection_TENDMLS(CPoint3F *pos, CVector3F *in, float *lastDTI, CVector3F *out);
	void GetFilteredDTI_TENDMLS(CPoint3F *pos, float *inDTI, float *outDTI);


private:
	void FiberTracking_FACT(const char *fiberPathname, const uint scale);
	void FiberTracking_STTMLS(const char *fiberPathname, const uint scale);

private:
	void ComputeFA();
	void CreateMLSSamplePoints();

	void DTIFilter_MLS(CPoint3F *pos, float *out_DTI, float *in_DTI);
	void DTIFilter_TENDMLS(CPoint3F *pos, float *out_DTI, float *in_DTI);
	void TrackingOneFiber_FACT(Seed *seed, CFiber *fiber);

private:
	void SaveFibers(const CFiber *fibers, const uint w, const uint h, 
		const char *confPathname, const char *vertPathname);

private:
	int ComputeIndex(CPoint3F *p) {
		uint xx = uint(p->m_x + 0.5f);
		uint yy = uint(p->m_y + 0.5f);
		uint zz = uint(p->m_z + 0.5f);

		int index = (int)((zz * m_height + yy) * m_width + xx);
		return index;
	}
};

#endif
